function ivp3

%  solve thermokinetic IVP using MATLAB routines

clear *
clf

ep=0.1;
tmax=8;

[t,y] = ode45(@rhs,[0 tmax],[0 0]);

% get(gcf)
set(gcf,'Position', [1281 735 577 197]);

tt=linspace(0,tmax,800);
for it=1:800
	u0(it)=1-exp(-tt(it));
	u1(it)=u0(it)+ep*( 0.5*(tt(it)^2+2*tt(it)-4)*exp(-tt(it)) + (2+tt(it))*exp(-2*tt(it))) ;
	q0(it)=1-(1+tt(it))*exp(-tt(it));
	q1(it)=q0(it)+ep*( (tt(it)^3-18*tt(it)+30)*exp(-tt(it))/6 - (2*tt(it)+5)*exp(-2*tt(it))) ;
end;

hold on
box on
yy=y(:,1);
plot(t,yy)
plot(tt,u0,'--r','LineWidth',1)
plot(tt,u1,'-.g','LineWidth',1)
%axis([0 2.5 0 0.601])
%title('\epsilon = 0.01','FontSize',14,'FontWeight','bold')
%legend(' Exact',' Asymptotic (1 term)',' Asymptotic (2 term)','Location','South');
legend(' Exact',' 1 term',' 2 terms',4);
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  
hold off

figure
hold on
set(gcf,'Position', [1283 465 577 197]);
box on
yy=y(:,2);
plot(t,yy)
plot(tt,q0,'--r','LineWidth',1)
plot(tt,q1,'-.k','LineWidth',1)
axis([0 8 0 1.08])
legend(' Numerical',' 1 term',' 2 terms',4);
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  
hold off

function dy=rhs(t,y)
dy=zeros(2,1);
ep=0.1;
dy(1)=1-y(1)*exp(ep*(y(2)-1));
dy(2)=y(1)*exp(ep*(y(2)-1)) - y(2);











